function [elements,processInfo] = CalculateMaterialResponse...
    (elements,model,processInfo)
ne = model.countElements();
dim = 2;
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    nn = f.countModes();
    ng = e.getNumGaussPoint();
    mE = elements(i).E;
    mNu = elements(i).nu;
    for j=1:ng
        xi = e.getGaussPoint(j-1);
        %compute strain
        B = zeros(3,dim*nn);
        B = e.computeB(f,B,xi);
        epsilon = B*ue;
        %compute trial deviatoric strain tensor
        vol_strain = epsilon(1) + epsilon(2);
        etrial(1) = epsilon(1) - vol_strain/3;
        etrial(2) = epsilon(2) - vol_strain/3;
        etrial(3) = epsilon(3); %engineering notation
        etrial(4) = -vol_strain/3;
        %compute trial elastic strain
        eedtrial = (etrial' - elements(i).ep_old(:,j));
        %compute trial deviatoric stress
        mu = 0.5*mE/(1+mNu);
        strial(1) = 2*mu*eedtrial(1);
        strial(2) = 2*mu*eedtrial(2);
        strial(3) = mu*eedtrial(3);
        strial(4) = 2*mu*eedtrial(4);
        %compute trial relative stress
        xitrial = strial' - elements(i).beta_old(:,j);
        normxi = sqrt(xitrial(1)^2 + xitrial(2)^2 +...
            2*xitrial(3)^2 + xitrial(4)^2);
        %check yield condition
        yield = computeIsotropicHardeningModulus(...
            elements(i),elements(i).alpha_old(j));
        ftrial = normxi - sqrt(2/3)*yield;
        if ftrial<0
            %update internal variables as trial state
            elements(i).alpha(j) = elements(i).alpha_old(j);
            elements(i).beta(:,j) = elements(i).beta_old(:,j);
            elements(i).ep(:,j) = elements(i).ep_old(:,j);
            %update constitutive matrix
            C = C_plane_strain(mE,mNu);
            elements(i).constitutiveMatrix(:,:,j) = C;
            %compute trial elastic strain
            eetrial = epsilon - elements(i).ep(1:3,j);
            %compute stress
            elements(i).stress(:,j) = C*eetrial;
        else
            %compute n
            n = xitrial/normxi;
            %compute delta gamma
            deltaGamma = computeDeltaLambda(elements(i),normxi,...
                elements(i).alpha_old(j));
            %compute alpha
            elements(i).alpha(j) = elements(i).alpha_old(j)...
                + sqrt(2/3)*deltaGamma;
            %compute beta
            H1 = computeKinematicHardeningModulus(...
                elements(i),elements(i).alpha_old(j));
            H2 = computeKinematicHardeningModulus(...
                elements(i),elements(i).alpha(j));
            elements(i).beta(:,j) = elements(i).beta_old(:,j)...
                + sqrt(2/3)*(H2-H1)*n;
            %compute plastic strain
            elements(i).ep(:,j) = elements(i).ep_old(:,j) +...
                deltaGamma*n;
            %compute stress
            kappa = mE/(3*(1-2*mNu));
            sigma(1) = kappa*vol_strain+strial(1)-...
                2*mu*deltaGamma*n(1);
            sigma(2) = kappa*vol_strain+strial(2)-...
                2*mu*deltaGamma*n(2);
            sigma(3) = strial(3) - 2*mu*deltaGamma*n(3);
            elements(i).stress(:,j) = sigma';
            %compute elastoplastic tangent
            theta = 1-2*mu*deltaGamma/normxi;
            dH = computeKinematicHardeningSlope(...
                elements(i),elements(i).alpha(j));
            dK = computeIsotropicHardeningSlope(...
                elements(i),elements(i).alpha(j));
            thetabar = 1/(1+(dH+dK)/(3*mu)) - (1-theta);
            onebyone = [1 1 0;1 1 0;0 0 0];
            rank4identity = [1 0 0;0 1 0;0 0 1/2];
            nbyn = [n(1)*n(1) n(1)*n(2) n(1)*n(3);
                n(2)*n(1) n(2)*n(2) n(2)*n(3);
                n(3)*n(1) n(3)*n(2) n(3)*n(3)];
            elements(i).constitutiveMatrix(:,:,j) = ...
                kappa*onebyone + 2*mu*theta*(rank4identity-...
                1/3*onebyone) - 2*mu*thetabar*nbyn;
        end
        
    end
end
end
